clear; clc; close all;

k = 3;
h = 0.01;
T = 1000;            % 积分总时间，较长保证轨迹稳定
transient = 600;     % 丢弃暂态
N = round(T / h);
N_trans = round(transient / h);

% m, n 取值集合
m_list = [-1, 0, 1];
n_list = [-1, 0, 1];

% 颜色设置，对应9组轨迹
colors = [
    0, 0, 0;          % 黑色 (0,0)
    0, 0, 1;          % 蓝色  (0,1)
    0, 1, 0;          % 绿色  (0,-1)
    0.5, 0, 0.5;      % 紫色  (1,1)
    1, 0.5, 0;        % 橙色  (1,0)
    1, 0, 1;          % 品红  (1,-1)
    1, 0, 0;          % 红色  (-1,1)
    0, 1, 1;          % 青色  (-1,0)
    0.6, 0.3, 0       % 棕色  (-1,-1)
];

figure; hold on; box on;

idx_color = 1;
for m = m_list
    for n = n_list
        mu = pi + 2 * m * pi;
        eta = 1 + 2 * n * pi;

        % 初始条件
        x = [1e-9; 0; 0; mu; eta];
        f = @(x) system_ode(x, k);

        % 丢弃暂态
        for t = 1:N_trans
            x = RK4_step(f, x, h);
        end

        % 采集轨迹点
        traj_u = zeros(N - N_trans, 1);
        traj_v = zeros(N - N_trans, 1);
        for t = 1:(N - N_trans)
            x = RK4_step(f, x, h);
            traj_u(t) = x(4);
            traj_v(t) = x(5);
        end

        % 绘制轨迹点
        plot(traj_u/pi, traj_v/pi, '.', 'Color', colors(idx_color,:), 'MarkerSize', 1);

        idx_color = idx_color + 1;
    end
end

xlabel('u', 'FontSize', 12);
ylabel('v', 'FontSize', 12);
grid on;

xticks([-2 -1 0 1 2 3 4]);
xticklabels({'-2\pi', '-\pi', '0', '\pi', '2\pi', '3\pi', '4\pi'});
yticks([-3 -2 -1 0 1 2 3]);
yticklabels({'-3\pi', '-2\pi', '-\pi', '0', '\pi', '2\pi', '3\pi'});

function dx = system_ode(x, k)
    % 状态变量 x = [x; y; z; u; v]
    dx = zeros(5,1);

    dx(1) = x(2) + x(3) - k * cos(x(5)) * x(2); % \dot{x}
    dx(2) = -x(1) + x(3);                        % \dot{y}
    dx(3) = -x(1) - x(3) + k * cos(x(4)) * x(1); % \dot{z}
    dx(4) = x(1);                               % \dot{u}
    dx(5) = x(2);                               % \dot{v}
end

%% 子函数：RK4一步积分
function x_next = RK4_step(f, x, h)
    k1 = f(x);
    k2 = f(x + h/2 * k1);
    k3 = f(x + h/2 * k2);
    k4 = f(x + h * k3);
    x_next = x + h/6*(k1 + 2*k2 + 2*k3 + k4);
end
